Process
- make buffer (400 m radius circles around well, which covers about 50 ha of land) for each well
source("./GitControlled/Codes/0_1_ls_packages.R")
source("./GitControlled/Codes/0_0_functions.R")
# /*===== point data for well location =====*/
data_w_LR_TB <- readRDS("./Shared/Data/WaterAnalysis/data_w_LR_TB.rds")
well_sf <- data_w_LR_TB %>%
unique(.,by="wellid") %>%
.[,.(wellid, longdd, latdd)] %>%
# --- the geographic coordinate in the well data is NAD83 (espg=4269)--- #
st_as_sf(., coords = c("longdd","latdd"), crs = 4269)
# /*===== make buffer around each well =====*/
well_buffer_sf <- well_sf %>%
#--- project to WGS UTM 14 (Cartesian 2D CS covering NE)---#
st_transform(32614) %>%
# --- make buffers--- #
# 400 m radius circles around well, which covers about 50 ha of land
st_buffer(., dist = 400)
ggplot()+
geom_sf(data=well_sf, size=0.5) +
geom_sf(data=well_buffer_sf, color="red", fill=NA)

well_buffer_sp <- as(well_buffer_sf, "Spatial")
- Then, apply
get_ssurgo_props() for each well buffer.
get_ssurgo_props <- function(field, vars, summarize = FALSE) {
# field=unique_well_sp; vars = c("sandtotal_r", "claytotal_r")
# Get SSURGO mukeys for polygon intersection
ssurgo_geom <-
SDA_spatialQuery(
geom = field,
what = "geom",
db = "SSURGO",
geomIntersection = TRUE
) %>%
st_as_sf()
mutate(
area = as.numeric(st_area(.)),
area_weight = area / sum(area)
)
# Get soil properties for each mukey
mukeydata <-
get_SDA_property(
property = vars,
method = "Weighted Average",
mukeys = ssurgo_geom$mukey,
top_depth = 0,
bottom_depth = 150
)
ssurgo_data <- left_join(ssurgo_geom, mukeydata, by = "mukey")
if (summarize == TRUE) {
ssurgo_data_sum <-
ssurgo_data %>%
data.table() %>%
.[,
lapply(.SD, weighted.mean, w = area_weight),
.SDcols = vars
]
return(ssurgo_data_sum)
} else {
return(ssurgo_data)
}
}
Example: SSURGO polygon data for an example field
library(soilDB)
library(aqp)
library(sp)
ex_well_bf <- well_buffer_sf[5,]
ex_well_bf_sp <- as(ex_well_bf, "Spatial")
ssurgo_geom <-
SDA_spatialQuery(
geom = ex_well_bf_sp,
what = "geom",
db = "SSURGO",
geomIntersection = TRUE
) %>%
st_as_sf()%>%
dplyr::mutate(
area = as.numeric(st_area(.)),
area_weight = area / sum(area)
)
# /*===== visualization =====*/
ggplot()+
geom_sf(data=ssurgo_geom, aes(fill=factor(mukey)), color="blue", alpha=0.5)+
geom_sf(data= well_sf[5,], fill=NA)
